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SYMBOLS 


A n load measurements in air, n = 1,2,3 

a model coefficient to be determined empirically, 

a = fl(a) 

Crf drag coefficient 

Cj lift coefficient 

C m quarter-chord pitching-moment coefficient 

c airfoil chord, m 

e model coefficient to be determined empirically, 

e = e(a) 

F aerodynamic load coefficient , F = F{r) 

F x component of F sufficient for linear range, F x = F l (r) 

F 2 supplement to F x required for nonlinear range, 
F 2 =F 2 (t) 

Fj linear extrapolation of the static load curve , Fj = F^a) 

F $ static load , F $ - F s (ol) 

F mean value in Taylor expansion for load, F = F(t) 

F modulus of first harmonic in Taylor expansion for 
load , F = F(t) 

fi forcing function on F x defined by second member in 
equation (2) 

f 2 forcing function on F 2 defined by second member in 
equation (3) 

k reduced frequency of oscillation, k = cjc/IU ^ 

Re Reynolds number, U^c/v 

r model coefficient to be determined empirically, 

r= r( a) 


s model coefficient to be determined empirically, 

s = s(a) 

t time, sec 

U ^ free-stream velocity, m/sec 

W n load measurements in water, n- 1,2,3 

a airfoil incidence, deg 

a 0 mean value of pitch oscillation, deg 

ax amplitude of pitch oscillation, deg 

A difference between linear and static load curves, Fj - F s 

e unity -step function 

X model coefficient to be determined empirically, 

X = X(a) 

v kinematic viscosity, m 2 /sec 

\ mn coefficients in mth cubic spline equation, n = 1,2,3 
p density, kg/m 3 

a model coefficient to be determined empirically, 

a = a(a) 

r reduced time, r = c ot/k 

co frequency of oscillation in pitch, rad/sec 

(*) derivative with respect to time, d/dt 

(**) second derivative with respect to time, 3 2 /dt 2 

4 imaginary part of quantity 

N total number of increments in one cycle 

<R real part of quantity 
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APPLICATION OF THE ONERA MODEL OF DYNAMIC STALL 


K. W. McAlister, 0. Lambert,* and D. Petotf 

Aeromechanics Laboratory, U.S. Army Research and Technology Laboratories, AVSCOM 


SUMMARY 


A semiempirical model , developed at the Office National D’Etudes et de Recherches Aerospatiales 
(ONERA), to predict the unsteady loads on an airfoil that is experiencing dynamic stall , is investigated . This 
study describes the math model from an engineering point of view , demonstrates the procedure for obtaining 
various empirical parameters, and compares the loads predicted by the model with those obtained in the 
experiment. The procedure is found to be straightforward , and the final calculations are observed to be in 
qualitative agreement with the experimental results. Comparisons between calculations and measurements 
also indicate that a decrease in accuracy results when the values of both the reduced frequency and the 
amplitude of oscillation are large. Potential quantitative improvements in the accuracy of the calculations 
are discussed in terms of accounting for both the hysteresis in the static data and the effects of stall delay in 
the governing equations. 


INTRODUCTION 


When a helicopter is in forward flight, the rotor blade 
must undergo a cyclic variation in incidence in order to 
balance the lift developed both over the advancing and 
retreating quadrants (to prevent roll), as well as over the 
forward and rearward quadrants (to prevent pitch). Increas- 
ing the flight speed of the helicopter increases the asymmetry 
in the local velocity distribution between the advancing and 
retreating sides. In order to eliminate the roll moment pro- 
duced by this asymmetry, the incidence of the rotor blade 
must be decreased on the advancing side and increased on the 
retreating side. Clearly there is a limit to how much the inci- 
dence can be increased without causing the flow to separate 
from the blade. Nevertheless, conditions often exist when 
this separation boundary is briefly penetrated, giving rise to a 
phenomenon called dynamic stall. 

The characteristics of dynamic stall are strongly influ- 
enced by the time-dependent nature of the viscous region 
surrounding the airfoil. Although flow reversal may have 
progressed over most of the upper surface of the airfoil 
during a rapid increase in incidence, the boundary layer will 
normally remain attached for angles of incidence well 
beyond the stall angle observed in a steady flow environment 
(refs. 1 and 2). The onset of stall is initiated by the growth 
and passage of a vortex over the airfoil. By the time the 
vortex has reached the trailing edge, the resulting unsteady 
flow values for the lift, drag, and pitching moment on the 
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airfoil may have doubled their maximum steady flow values 
(ref. 3). As the vortex is swept into the wake of the airfoil, a 
sudden reversal in the lift and drag loads occurs, thereby 
creating a significant source of vibration on the helicopter. 
Another damaging aspect of dynamic stall is the rapid growth 
and decline of a nose-down pitching moment. This results 
from a rearward movement of the center of pressure which 
accompanies the passage of the vortex over the airfoil and 
into the wake. This impulsive character of the pitching 
moment acts as a strong forcing function on the aeroelastic 
stability of the rotor blade (ref. 4). A potentially dangerous 
situation may therefore develop if the interaction between 
the blade and the surrounding air results in oscillations that 
are negatively damped. Since rotor blades are typically both 
slender and flexible, this unstable condition, known as “stall 
flutter,” can be especially threatening to the safety of the 
helicopter. 

In order to calculate the performance boundaries of a 
potential rotor design, it is essential that the mathematical 
formulation adequately models the effects of dynamic stall. 
Unfortunately, a closed-form solution does not exist. In fact, 
all of the currently available methods employ some form of 
approximation or empiricism, and are also normally 
restricted to two-dimensional flows (ref. 5). Many of these 
prediction techniques are based directly on a recognition of 
the global attributes contained in the unsteady force and 
moment responses that are observed in the angle-of-attack 
domain. As a consequence, expressions have been devised by 
numerous investigators that describe explicitly the “time- 
delay” character of the various events embodying dynamic 
stall. A particularly noteworthy representative of this 
approach is presented in reference 6. Two time constants are 
featured in this model; one describes the time delay after 



exceeding the static stall angle and before a vortex is shed 
from the leading edge of the airfoil, and the other describes 
the time required for this vortex to reach the trailing edge. 

Generally speaking, most prediction techniques for 
dynamic stall have been successful only within the limits of 
the data from which they were fabricated. A model that is 
very much dependent on its data base is described in refer- 
ence 7. However, a fairly extensive range of data was con- 
sidered initially, and allowance was made for expansion. 
This method is based on a set of algebraic equations contain- 
ing parameters that must be determined from the available 
synthesized experimental data. Published results show that 
the procedure accurately reconstructs the aerodynamic loads 
that occur during dynamic stall over a wide range of condi- 
tions. An important factor contributing to the success of 
this method lies in the particular set of dynamic parameters 
that were postulated. The accuracy of this approach is depen- 
dent on the correctness of three semi-empirical expressions 
that describe (a) the airfoil incidence when moment stall 
occurs, (b) the dimensionless time when the stall vortex 
reaches the trailing edge, and (c) the incidence when the flow 
reattaches to the surface of the airfoil. 

In contrast to the methods that attempt to duplicate the 
effects of dynamic stall, a unique method was developed at 
ONERA that utilizes the characteristics of differential equa- 
tions to directly simulate the aerodynamic responses in the 
time domain (ref. 8). Although other techniques may render 
more accurate predictions of rotor blade loading, the primary 
advantage of the ONERA model is that the governing system 
of equations can be readily linearized, therefore making it 
well suited for inclusion in analyses of rotor stability (refs. 9 
and 10). Certain aspects of the model are continuing to 
undergo refinement; however, the fundamental concept 
appears to be well established (refs. 1 1 and 12), and will not 
be restated in detail here. Instead, the scope of this presenta- 
tion will be to describe the model from an engineering point 
of view, to demonstrate the attainment of various empirical 
parameters, and to compare the loads predicted by the model 
with those obtained from an experiment. 


DESCRIPTION OF THE MATH MODEL 


General Equations 

A fundamental assumption is made that the aerodynamic 
loads can be determined from a set of nonlinear transfer 
functions containing input variables that describe the motion 
of the airfoil. Considering the operational and structural 
environments that are typical for helicopter rotor blades, it 
is further assumed that all input and output variables are 
first-order small quantities and that the coupling can be 
neglected between the chord force and either the normal 
force or the pitching moment. Perhaps the most restrictive 


assumption to be imposed, in light of the model’s applicabil- 
ity to separated flows, is that the instantaneous force and 
moment loads do not depart too greatly from their steady 
flow values. This condition is necessary so that the coeffi- 
cients appearing in the governing equations reduce to func- 
tions of only velocity and angle of attack for a given airfoil. 
Additional simplifying assumptions are then made concern- 
ing a reduction in the order of the equations and the elimina- 
tion of certain cross coupling coefficients. The validity of 
these assumptions could not be demonstrated analytically, 
and had to be confirmed by experiment. A significant out- 
come from this simplification is that the model can now be 
used to evaluate any given load, independent of the others. 

Evaluating the model at different amplitudes and frequen- 
cies provided the necessary guidance in reaching the final 
form of the equations. Airfoil oscillations below stall indi- 
cated that the loads could be represented by a first-order 
equation, having a real negative pole; whereas oscillations 
beyond stall produced loads that required a second-order 
equation representation, having two complex conjugate 
poles. In view of these observations, the single-equation 
formulation for each load was abandoned. Instead, each load 
would be divided into two components, one governed by the 
first-order equation and the other by the second-order equa- 
tion. Letting the function F denote the total aerodynamic 
load of interest, the governing equations become 

F = F l +F 2 (1) 

F x + XFi = XFf + (As + a)d + sot (2) 

F 2 + aF 2 + rF 2 = -(] rA + e A) (3) 

where the coefficients A, s , a, a, r, and e will be treated as 
functions of a only, the instantaneous incidence of the air- 
foil. Strictly speaking, these coefficients are also dependent 
on the free-stream velocity (or alternatively on the Mach 
number) and the profile of the airfoil. However, to illustrate 
the application of this technique, it is sufficient to consider 
an incompressible flow so that a becomes the only param- 
eter. Since these coefficients actually represent time deriva- 
tives (no longer explicitly apparent), and must be determined 
experimentally, they can be obtained by performing small- 
amplitude oscillations around discrete values of a. The range 
of angles over which the coefficients must be specified is 
dictated by the large-amplitude cases to be calculated by the 
model. The variables F j and A, also functions of a, are com- 
pletely determined by the static behavior of the airfoil. F j 
denotes a linear extrapolation of the static load curve and 
A is defined to be the difference between this extrapolation 
and the actual static curve. Using a hypothetical lift response 
as an example, the relation between various terms is illus- 
trated in figure 1. 
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As the frequency of the airfoil oscillation diminishes in 
the limit, all time derivatives vanish and equations (l)-(3) 
reduce to 

lim F = Fj - A = F 

T~* 0 S 

where F s is the static load response. If the airfoil motion is 
unsteady, but remains entirely within the linear range, then 
A = 0 and the load is completely determined by the solution 
of equation (2) for F x . 

Small-Amplitude Equations 


with a can be ascribed. In this case, the function A(a) can be 
written as 




Although the variation of F j with a is simply a straight line, 
regardless of the amplitude of oscillation, it will nevertheless 
be expressed in the same form as for A; so that 



To determine the relationships between the six coeffi- 
cients and a, an experiment is required in which the airfoil is 
made to undergo pitching oscillations at different frequencies 
and mean angles, and for values above and below stall. 
Although the pitching motion can be random (ref. 11), only 
harmonic variations in incidence will be considered. Letting 
a 0 denote the mean angle, ct l the amplitude, and k the 
reduced frequency of oscillation, the incidence can be 
generally expressed as 

ol = a 0 + (R [cti e^ r ] (4) 

and is not restricted to small-amplitude motion. However, if 
the amplitude of the oscillation is small, say ot x < 1°, then 
the corresponding load can be approximated by 

F = F + (R[F e ikT ] (5) 

where F represents the mean value of the load, not necessar- 
ily equal to F , and Fis a complex function representing the 
modulus of the first harmonic of the load. This is only an 
approximate expression for the load since the forcing func- 
tion in equation (3), and the coefficients in both equa- 
tions (2) and (3), are nonlinear functions of a. As such, 
equation (5) is more correctly written with higher order 
terms; however, experiments have shown that when the 
amplitude of the oscillation is sufficiently small, the load 
response is nearly elliptical so that the first-harmonic repre- 
sentation is acceptable . 

Given the experimental observation that for small- 
amplitude oscillations a first-order input in a produces a first- 
order output in F, the governing equations (l)-(3) can be 
simplified to a system with locally constant coefficients and 
with forcing functions that are first-order harmonic relative 
to a. The forcing function in equation (2) already has a first- 
harmonic form since it is linear in a and its derivatives. The 
same is not generally true of the forcing function in equa- 
tion (3) since the term A is normally quite nonlinear with a. 
However, given the restriction that the oscillations will be 
small in amplitude means that a piecewise -linear variation 


Reclaiming the imaginary parts of a and F in equa- 
tions (4) and (5) above, the dependent and independent 
variables become 

a = a 0 + a x e*^ T (8) 

Ft =Fi +Fi e ikT (9) 

F 2 =F 2 + F 2 e ikT (10) 

Substituting the quantities given in equations (6)-(10) into 
the governing equations (l)-(3) yields 

F = F, +F 2 +(F 1 +F 2 )e ikT (11) 

- .. r dF, 

XFj + (X + ik)F x e tkT = XF,| + X — 

a ° L da «o 

+ ik(Xs + a) - a, 4 kj (12) 

rf' 2 + (r + iak - k 2 )F 2 z‘ kr = -r A| - (r + iek)a l e‘ kr — 

° da 

(13) 

Equations (12) and (13) both state that two equalities must 
be satisfied; the terms of one equality are steady while those 
of the other are unsteady (identified by their product with 
the complex potential term). The equalities that are com- 
posed of steady terms are: 


and 

■ - a 'o 0 

Summing the above two expressions to obtain the mean 
value for the load gives 
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F = F } 


Pot 0 


. A| =F 
ol 0 s 


(14) 


and states that when the amplitude of the oscillation is small, 
the mean load will be the same as its static value. In actual- 
ity, the mean value will depart slightly from its static value as 
the frequency of oscillation is increased; however, the linear- 
ization of the equations prevents this behavior from being 
reproduced. Taking similar steps to arrive at the unsteady 
portion of the load results in 


F 


da 


+ ik(\s + a) - sk 




2 


(X + ik ) 


(r + iek) 


dA 

da 


'<*o 


(r + iak - k 2 ) 


(15) 


This relation represents the transfer function of the oscillat- 
ing portion of the load (F) relative to the input (c*i ). Separat- 
ing the real and imaginary components of equation (15) 
yields 


« 


F 

JF l 

Qi 

da 




(x 2 + r ) 



k 2 (r-ae)-r 2 dA 
(k 2 - r) 2 + (ak) 2 da 


la 0 


(16) 


and 
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F 


«i 


= ks + 


k\ 

(X 2 + k 2 ) 



+ ek(k 2 - /*) + akr dA 
+ {k 2 - r) 2 +(akf da 


(17) 


and for large values of k, assuming that k » X, equa- 
tion (16) shows that o is the asymptotic value of the real 
part of the load; similarly, equation (17) shows that s is the 
asymptotic value of the rate of change of the imaginary part 
of the load with respect to the reduced frequency. In other 
words, 


lim — = o+iks (18) 

k^oo a x 

The small-amplitude equations are now in a convenient 
form to evaluate the six coefficients, and they can be applied 
over the entire incidence range, both above and below stall. 
To evaluate these coefficients, an experiment must be per- 
formed to obtain the load measurements during small- 
amplitude oscillations. The measurements are then Fourier 


analyzed to determine the real and imaginary components of 
the first harmonic. The coefficients, which are taken to be 
constant during the oscillation around any given mean angle, 
must assume values as required to satisfy the equality 
between the real and imaginary load measurements and those 
described by equations (16) and (17). After obtaining these 
coefficients, and knowing the static behavior of the airfoil, 
the governing equations (l)-(3) can be solved to obtain the 
load. This entire sequence is diagrammed in figure 2. 


DESCRIPTION OF THE EXPERIMENT 


The experiment was conducted in the 4000-liter, closed- 
circuit facility at the Aeromechanics Laboratory Water 
Tunnel, Ames Research Center (fig. 3). The test section is 
21 cm wide, 31 cm high, and extends horizontally a distance 
of 86 cm. The airfoil selected for this study was a Boeing- 
Vertol VR-7, having a two-dimensional planform with a 
chord of 10 cm. The airfoil was positioned so that it spanned 
the width of the test section to within 0.015 cm of either 
side. In order to minimize the moment of inertia about its 
pitch axis, the airfoil was cast from a lightweight epoxy resin 
around a metal spar. The pitch axis was placed at the quarter- 
chord location. When installed, the spar of the airfoil 
extended through the test-section windows and was sup- 
ported by lift and drag transducers on both sides (fig. 4). One 
end of the spar was adjoined to an instrumented drive shaft 
through a torsionally stiff coupling so that airfoil incidence 
and pitching moment could be measured. Frictional 
moments imparted by the support bearings and seals were 
also measured and later treated as dynamic-load tares. 

Electrical instrumentation consisted of transducers for 
the measurement of airfoil incidence, lift (both sides), drag 
(both sides), total pitching moment, and the bearing and 
seal moments (both sides). After amplification, these signals 
were either appropriately summed (i.e., total pitching 
moment minus both frictional moments) and displayed on 
local monitors; or they were transmitted to a remote data 
acquisition system where they were digitized, averaged, and 
stored for later processing. Digitizing and ensemble averaging 
was based on two additional signals: a 360/rev pulse train 
that was synchronous with cot, and a 1/rev pulse that was 
synchronous with the beginning of each cycle of airfoil 
oscillation. On-line monitoring for smoothness of the ensem- 
ble average of a particular load provided the basis to termi- 
nate the data acquisition; and the number of cycles used in 
generating the average, therefore, was dependent on the 
extent of the nonperiodic content of the signal. It is esti- 
mated that the incidence of the airfoil could be measured to 
an accuracy of 0.2°. Lift and drag measurements are con- 
sidered to be accurate to 0.01 N and the pitching moments 
to 0.002 Nun. 
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With the airfoil set at zero incidence, the speed of the 
water was adjusted to produce a dynamic pressure of 
689 N/m 2 (0.10 lb //in. 2 ). This pressure corresponds to a 
Reynolds number of Re = 120,000 based on the chord of the 
airfoil. The tunnel was then operated at a fixed drive speed 
for the duration of the experiment. Some reduction in tunnel 
speed is thought to have occurred when the airfoil was 
stalled; however, no attempt was made to either measure or 
account for this degradation. 

In order to obtain the asymptotic character of the real 
and imaginary components of the load at high frequency, say 
at k = 1.2, an oscillation frequency of around 4.5 Hz would 
be required. Inertial loads might, therefore, become a prob- 
lem. The inertial loads are potentially less serious in water 
than in air because of the relatively low density ratio 

Pmodel/Pwater compared with P m odel/Pair- However, at the 
beginning of this test it was not known whether the inertial 
effects on the hydrodynamic load responses could be 
neglected. Furthermore, there was also concern about 
gravitational effects (buoyancy) and the possibility of 
variable support loads (due to misalignments) during angle 
changes. The contribution of all such loads to the balance 
measurements is summarized schematically in figure 5 . 

To obtain the desired hydrodynamic load, the gravita- 
tional, support, and inertial loads must be removed from the 
unsteady load measurements. This can be accomplished 
during the data reduction phase, provided that the appro- 
priate quantities are measured at the time of the experiment. 
The inertial effects can be accounted for by lowering the 
water in the test section and performing the same unsteady 
load measurements in air. The gravitational and support loads 
can be accounted for by performing quasi-static measure- 
ments in water (at zero flow) as well as in air. Having 
recorded these different results, the final hydrodynamic load 
can be calculated by following the procedure outlined in 
figure 6. 

DISCUSSION OF RESULTS 


Load Measurements 

The variables Fjfot) and A(a) are both essential elements in 
the math model, and they can easily be determined from the 
static behavior of the airfoil. Their variation with a must be 
established over the entire incidence range for which large- 
amplitude responses are to be calculated. Static measure- 
ments of lift, drag, and pitching moment are shown in 
figure 7. In order to evaluate the capability of the model for 
calculating loads under dynamic stall conditions, measure- 
ments were also made with the airfoil undergoing large - 
amplitude oscillations over a range of reduced frequencies 
(figs. 8-14). 


To determine the six coefficients of the math model, 
small-amplitude oscillations in pitch are required around a 
selection of mean angles in the linear and nonlinear range. 
For some airfoils, the static stall angle marks the division 
between these two ranges; however, as can be seen from 
figure 7, the nonlinear range in the present case begins well 
below the stall angle. Within the linear range, three of 
the coefficients become inconsequential because A = 0; 
while among the remaining three coefficients to be evaluated, 
two are usually found to be constant (ref. 12). As a result 
of this reduction in complexity, data are needed for only 
a few mean angles in the linear range. In the present case 
the linear range appears to be so brief that little advantage 
can be realized in the number of mean angles to be 
required. 

In keeping with the assumptions under which equa- 
tions (16) and (17) were derived, the amplitude of oscillation 
was restricted to = 0.5°. Measurements were made at 
13 mean angles (a 0 = 0°, 1.5°, 3°, 5°, 7°, 9°, 1 1°, 13°, 15°, 
17°, 19°, 21°, and 23°) and 11 reduced frequencies 
(k = 0.025, 0.05, 0.10, 0.15, 0.2, 0.3, 0.4, 0.5, 0.7, 0.9, 
and 1.2). Excluding the lowest-frequency case, these small- 
amplitude results are presented in figures 15(a)-15(m) for 
the lift, drag, and pitching moment. Also included are the 
first-harmonic equivalents of these data (shown as dashed 
lines). The high frequency component of the lift and 
pitching-moment signals, especially apparent for values of 
incidence below 1.5°, is due to Karman-vortex shedding from 
the airfoil. The frequency of this vortex shedding appears to 
be around 47 Hz. The first-harmonic curves appear to com- 
pare well with the actual data, and therefore will be accepted 
in the following analysis as an accurate replacement for the 
small-amplitude measurements. Examining these first- 
harmonic curves in more detail, the real and imaginary com- 
ponents of each load modulus are separated and shown in 
figure 16. For any given load, these components correspond 
to <R[F] and .?[F] in equations (16) and (17), and as such 
provide the basis on which the six coefficients can be 
determined. 


Coefficient Evaluation 

The most important step in finalizing the governing differ- 
ential equations for the loads on a given airfoil is the deter- 
mination of the equation parameters. Specification of the 
parameters F/(a) and A(a) evolves directly from the static 
load measurements, which are shown in figures 17 and 18. 
To smooth out any irregularities in the experimental data, 
and to provide a more well-behaved variation with incidence 
so that reasonably bounded derivatives could be calculated, 
the static curves have been fitted with piecewise -cubic 
splines. Both A and dA/da are derived using these curve fits. 

Evaluation of the six coefficients is less straightforward. 
Generally speaking, coefficients cannot be found which 
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( 20 ) 


precisely satisfy equations (16) and (17) over the frequency 
range desired, and at the same time culminate in a smooth 
dependence on incidence. Instead, solutions must be sought 
at each mean angle which satisfy the equations for (R[F] and 
#[F] in a least-squares sense. Furthermore, experience has 
shown that a particular order should be followed in establish- 
ing the set of coefficients. Since the primary intent of this 
presentation is to provide an example of the application of 
the ONERA model, the following discussion will be limited 
to the lift coefficient so that the method can be discussed in 
more detail and with greater clarity. 

Before evoking this procedure, it is useful to identify 
those coefficients that, from experience, have been found to 
be independent of incidence and whose values are essentially 
declarable by inspection. Considering the linear range, equa- 
tion (16) states that 


fl 


F_ 

0L\ 


dFj 
\ 2 — 
da 


+ ok 2 


<*o 
X 2 +* 2 


(19) 


This expression can be used to examine the bounds on 
<R[F/ai ] as A; is varied. On one extreme, k = 0 requires that 
= ( dFi/da ) ao . On the other extreme, letting 
k large suggests that (k[F/ai ] -*o. Although o is in general 
a function of a, its value can be readily approximated in the 
linear range. It can also be shown that the average value for 
<R[F/ai] over this extreme range of reduced frequencies is 
[(dFi/da) a + a]/2, and that this value is obtained when 
k = X. Accordingly, a second coefficient can be approximated 
in the linear range. Experience also indicates that the magni- 
tude of X has an effect on the manner in which the real and 
imaginary asymptotes are approached in the linear range; 
small values accelerate the approach and large values cause 
the approach to be more gradual (ref. 13). Additionally, 
some observations have been made about X in the nonlinear 
range. Experience has also shown that for low and medium 
reduced frequencies, F 2 is dominant over . Recalling from 
equation (15) that the first term of the second member 
represents F x and that the second term represents F 2 , it is 
evident from the composition of these two quantities that 
the value of X is of no consequence in this case. Furthermore, 
an examination of equation (18), which is valid at high 
reduced frequencies, shows that X is completely absent from 
the expression for F (which equals the sum of F Y and F 2 ). 
This implies that the value of X has a negligible effect on the 
math model for all frequencies in the nonlinear range, and 
that its value in the linear range can simply be extrapolated 
over the entire incidence range. As a result, it has been the 
practice at ONERA to simply accept a constant value for 
X. In keeping with this attitude, X can be approximated from 
the data for the real component of the lift at a 0 = 0° shown 
in figure 16(a): 


X = 0.25 (Va) 

Another coefficient that appears to vary little with inci- 
dence is 5 . Representing the asymptotic slope of the imagi- 
nary part of the load at high reduced frequencies, the experi- 
mental data presented in figure 16 confirm that s can indeed 
be approximated by a constant over nearly the entire inci- 
dence range. For those cases in the nonlinear range where the 
asymptotic value has not yet been reached at the highest 
reduced frequency under consideration (i.e., at a 0 = 23 ), it 
is common practice to accept the extrapolated value for s 
from the linear range. Taking an average of these measure- 
ments yields 


s - 0.12 (Va) (21) 

Determination of the remaining four coefficients is more 
laborious since they must be considered simultaneously, they 
typically all vary with incidence in some nonlinear fashion, 
and yet nowhere do they satisfy the governing equations 
precisely. The optimizing algorithm selected for this task is 
based on a numerical scheme that finds solutions to equa- 
tions (16) and (17) so as to produce the least accumulated 
disagreement with the measured load responses. Considering 
that these coefficients will vary mostly in the nonlinear 
range, and since A is a more direct measure of the departure 
from the linear curve than is a , all four coefficients will be 
considered as explicit functions of A instead of a. This 
implies that their values in the linear range (where a may 
vary greatly) will be invariant and equal to their values at 
A = 0. 

Exercising the optimization routine with all four coeffi- 
cients initially unknown yields an expected dispersion of 
solutions (fig. 19). The customary step at this point in most 
parameter-identification procedures is to recognize which 
feature or characteristic of the data is most prominent, and 
on which coefficient is this feature most dependent. Sensi- 
tivity tests have shown that the location of the peak values 
of both the real and imaginary responses is such a feature, 
and that it is strongly (but not solely) dependent on r. The 
optimization for r will therefore be examined first, and a 
curve defined that best describes this coefficient’s variation 
with A. Although a simple parabolic relationship might have 
sufficed in this case, the distribution of points for the other 
three coefficients indicated that a more flexible curve fit 
would be necessary. To be consistent, a piecewise -cubic 
spline was chosen for all curve fits. A curve is first faired 
through the data and then a set of manufactured data points 
is defined along this curve. On any given interval between 
these manufactured points, the curve is described by: 

ft" ~ + ~ a ) + ~ 3 

( 22 ) 
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where m represents the interval for which a m < a < a m+1 . 
Obviously, this curve fit is not unique, and a great deal of 
subjectivity can influence the weighting of various data 
points. In such cases, experience can become an important 
factor in the shaping of curves through the data. In any case, 
it may still be necessary to make certain adjustments to the 
coefficients after reviewing how the model calculates a large- 
amplitude case. For example, assume that several calculations 
are performed for an airfoil undergoing large-amplitude oscil- 
lations through stall. Furthermore, assume that the airfoil 
then returns to attached flow conditions with suspiciously 
little change in the hysteresis, in each case, as the reduced 
frequency is increased. This may suggest that the values of 
a and r in the linear range (hence their value at A = 0) need 
to be modified to allow less damping and more response in 
the equation for F 2 . This behavior could not have been fore- 
seen from small-amplitude tests, and perhaps not even from 
large-amplitude tests where the reduced frequency is small. 
However, for large-amplitude oscillations at moderate 
reduced frequencies, the time-history effects become suffi- 
ciently significant and can influence the loads even though 
the forcing function in the equation for F 2 may be locally 
small or zero. 

Although the curve for r can be regarded as provisional, 
and eventually may be modified, it will be considered as 
fixed for the present. With X, s , and r specified, the optimiza- 
tion process is repeated; but now only a, a , and e are con- 
sidered to be unknown. As can be observed in figure 20, the 
dispersion of points for the remaining three coefficients is 
considerably reduced, thereby the plotting of succeeding 
curve fits can be done with greater confidence. Representing 
the asymptotic value of the real part of the load, the data for 
a will be considered as the most reliable even though the 
maximum reduced frequency at each a 0 generally appeared 
to be too low to reveal the actual asymptote. A simple curve 
fit through the data for o would be especially infeasible in 
this case. For the special case when A = 0, the value for a 
used to evaluate X (recall when a = 0°) also applies. For 
values of A near zero, o is perhaps more dependent on a 
(explicitly) than on A; whereas for larger values of A, a very 
demonstrative dependence on A is observed. In fairing a 
curve through the data for o, those points for which A « 0 
(and which may eventually require some explicit dependence 
on a) were ignored and an expression similar to (22) was 
established for a(A). 

With only two unknowns remaining, a and e , the optimi- 
zation procedure yields a set of new values for each coeffi- 
cient with even less dispersion (fig. 21). Like a, both a and e 
exhibit a behavior around A « 0 that may require an explicit 
dependence on a if strict adherence is to be retained. How- 
ever, in keeping with the decision made for a, those points 
will be ignored for the present. Instead, the value of each 
coefficient at A = 0 will be that obtained by extending the 
curve from higher values of A. The coefficient a affects both 
the amplitude and width of the mid-frequency wave that is 


characteristic of the real and imaginary responses in the non- 
linear range (a 0 > 9°). The data in this range appear to be 
reasonably well behaved so that a curve similar to the one 
described by (22) can be established for <z(A). This curve, 
when extrapolated back to A = 0, yields a value for a that is 
substantially higher than what the optimizations (as well as 
the experimental observations) require. During a calculation 
for an airfoil at low incidence recovering from stall, this 
region will contribute significantly to the damping in the 
equation for F 2 . Therefore, some modification to the curve 
for a as A -» 0 can be anticipated. 

The optimization process can now be narrowed down to 
the evaluation of a single unknown, e, and the resulting dis- 
tribution of computed values for this coefficient is shown in 
figure 22. Although the variation of e with A is quite differ- 
ent after stall (A > 0.9) compared with its behavior over the 
remainder of the nonlinear range, the fitting of a piecewise- 
cubic spline to the data is straightforward. As can be seen 
from the governing equations (l)-(3), the coefficient e con- 
tributes only to the forcing function for F 2 . Aside from its 
impact on the time history of the solution as A -> 0, the fact 
that e appears as a product of A means that its value in the 
linear range is not important. The unknown e affects the 
amplitude of both the real and imaginary parts of the 
response in the nonlinear range as well as the phase of the 
response relative to the motion of the airfoil. Negative values 
of e cause a phase lag with respect to a; and as the magnitude 
of e increases, the phase lag tends toward 90°. The opposite 
trend occurs for positive values of e. 

To conclude the optimization process, the expressions 
obtained for the six coefficients will be coupled with equa- 
tions (16) and (17) to examine the real and imaginary com- 
ponents presently represented by the model. The focus of 
the examination will be on how well the model reproduces 
the real and imaginary components of the lift coefficient 
obtained in the experiment. Figures 23 and 24 show the 
comparison between the experimental values (discrete sym- 
bols) and the calculated values (solid curves) at the mean 
angles for which data exist. Excellent agreement is obtained 
through a 0 = 15°; however, the calculations resemble only 
the trend of the experiment for mean angles between 
17° < a 0 < 23°. That such good agreement exists between 
the calculations and the experiment for low values of a 0 (in 
spite of a disregard for many of the computed optimizations 
in the range A ^ 0) suggests that the small amplitude tests 
provide little information about the coefficients r, a, a, and e 
in this domain. The less satisfactory agreement between the 
calculations and the experiment for the higher values of a 0 
(A > 0.2) may be either due to an inadequacy in the model’s 
equations to account for a fully separated flow or due to the 
optimization process which assigned equal weighting to each 
of the experimental observations. In any case, a judgment on 
the consequence of this shortcoming will have to await the 
calculation of a large-amplitude case. The topographies of 
the real and imaginary components over the k-a 0 plane are 
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shown in figure 25 for the experimental data and in figure 26 
for the math model using the coefficient optimizations 
(recall fig. 22). 


Large-Amplitude Calculations 


Using the expressions for the six coefficients obtained 
above, the unsteady load for a prescribed airfoil motion can 
be calculated from equations (l)-(3). First, the governing 
system must be reduced to a set of first-order equations in 
order to use a standard differential-equation solver. The 
governing equations are readily transformed to the following: 


F=F i +F 2 

F x = ~XF x + \Fj + (Xs + o)ol + sa 


F 7 = F : 


> 


(23) 


F 3 = -aF 3 - rF 2 - (rA + e A) 


J 


Given a set of initial conditions, say Ffir = 0) = 0.0 
(i = 1, 2, 3), a time marching solution can be generated with 
dimensionless time r, the independent variable, expressed as 


(24) 

where N denotes the total number of increments in one 
cycle. The number of time steps required to reach a steady- 
state solution depends on the magnitude of the unsteady air- 
foil motion. Therefore, the higher the reduced frequency, A:, 
the greater the number of time steps required. 

By requiring that r = 0 when a = 0°, the sinusoidal motion 
of the airfoil can be described by 


a = a 0 - a x cos rk (25) 

The calculations are based on a mean angle and amplitude of 
oscillation of a 0 = Qa - 10° and a range of reduced frequen- 
cies from 0.002 < k < 0.25. These conditions were chosen 
because they challenge the math model to accurately predict 
the dynamic loads occurring during deep stall and with vary- 
ing degrees of overshoot (values above static Cj-max) and 
hysteresis. The results for the lift coefficient are shown in 
figures 27(a)-27(g), and a steady-state solution was reached 
during the first cycle of oscillation for all of the cases 
calculated. 

Each figure is subdivided into three subplots. The right- 
hand subplot displays the calculated components of the lift 
coefficient, F x and F 2 , as well as the measured static 
behavior of the airfoil (dashed line). The lower-left subplot 
shows the forcing functions f x and f 2 that appear in the 


computation of equations (2) and (3), respectively. The 
upper-left subplot shows the calculated value for the lift 
coefficient, F x + F 2 , along with the actual measured 
response (dashed line). The results given in these figures indi- 
cate that the math model does reproduce qualitatively the 
increases in both overshoot and hysteresis as the reduced 
frequency is increased. Even the slight surge in the lift 
response just prior to stall (as well as just prior to reattach- 
ment) is correctly predicted by the math model (i.e., see 
fig. 27(c)). This characteristic is believed to be attributable 
to the term A that appears in the forcing function for F 2 , 
and reflects the abrupt static stall that was actually measured 
in the experiment (recall figs. 17 and 18). However, the 
extent of both the overshoot and the hysteresis appears to be 
quantitatively incorrect. In addition, the calculated lift coef- 
ficient during the deep-stall phase of the motion (a < 0) is 
generally too low, and becomes even more so as the reduced 
frequency is increased. It appears that a decrease in accuracy 
results when the values of both the reduced frequency and 
the amplitude of oscillation are large. 

To overcome the quantitative disagreement in the lift 
overshoot between the calculations and the measurements, it 
appears that the math model must be amended to better 
account for the effects of stall delay. This eventuality was 
envisaged early on by the architects of this model (ref. 8). 
Consequently, an allowance for the observed delay in stall 
was provided for through the product of a unity-step func- 
tion, e, with the forcing function on F 2 in equation (3). The 
intent of this procedure is to disable the forcing function on 
F 2 during a certain period of time, dr, beginning when the 
static stall angle is exceeded. This dimensionless-time interval 
is defined by 



and typically has been set at about 5 r ^ 10 in recent applica- 
tions of the model (ref. 12). After satisfying this delay in 
time, the forcing function then assumes its current value. The 
time-delay concept (a prominent feature in ref. 6) is physi- 
cally based on the overshoot caused by the stall vortex. This 
vortex induces additional lift on the airfoil when passing over 
the airfoil. 

Since vortex shedding just prior to reattachment is gener- 
ally not obvious from available experimental data, it is 
unlikely that a similar time-delay factor applies to the stall 
recovery phase of the cycle. However, since a hysteresis 
clearly exists in the static response, it does seem appropriate 
to distinguish between a > 0 and a < 0 in terms of the defi- 
nitions for A (and therefore A). Doing so would automati- 
cally establish a persistence in the forcing function on F 2 and 
cause an extension of the hysteresis. A proper accounting for 
the stall-delay effect during a > 0, as well as the recognition 
of a separate A during a < 0, both have a potential for 
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improving substantially the quantitative accuracy of the 
math model. As such, these considerations need to be evalu- 
ated before proceeding to the modeling of the drag and 
moment responses. 

CONCLUSIONS 

A semiempirical model, developed at the Office National 
D’Etudes et de Recherches Aerospatiales (ONERA) to 
predict the unsteady loads on an airfoil that is experiencing 
dynamic stall, has been described. Calculations were per- 
formed for comparing with results from an experiment in the 
Aeromechanics Laboratory Water Tunnel. The static and 
dynamic data obtained were for a Boeing-Vertol VR-7 airfoil 
at a Reynolds number of Re = 120,000. Although lift, drag, 
and moment responses (both small and large amplitude) have 
been included in this study, the comparisons that have been 
discussed, including the following conclusions, are confined 
to the lift coefficient. 

1. The procedure, established at ONERA to identify the 
six coefficients in the math model, was found to be straight- 
forward and the trends consistent with earlier investigations 
when used with small-amplitude experimental data. 

2. Coefficient optimizations were found that produced 
excellent agreement between calculated and measured small- 
amplitude responses, as long as the mean angles were within 
0° < a 0 < 15°. For the higher mean angles that correspond 
to deep stall, 17° <a 0 <23°, the calculated small-amplitude 
responses based on optimized coefficients in this range were 
not in good agreement with measured values. 

3. The shape of the static-lift curve in the present experi- 
ment differs from the idealized version (as well as from mea- 
surements on other airfoils at higher Reynolds numbers) 
around which the math model was originally formulated. 


In the present case, the static curve becomes nonlinear prior 
to stall and deviates greatly from the classic tangent to the 
data near a = 0°. Furthermore, when stall does occur, it is 
quite abrupt. The extent to which this behavior affects the 
applicability of the model has not yet been resolved. 

4. A comparison between calculations and measurements 
for large -amplitude cases at various reduced frequencies 
shows that the math model does reproduce qualitatively the 
increases in both dynamic-lift overshoot (values above static 
Ci-max) and hysteresis as the frequencies are increased. Even 
the slight surge in the lift response just prior to stall (as well 
as just prior to reattachment) is correctly predicted by the 
math model. However, the extent of both the overshoot and 
the hysteresis appears to be quantitatively incorrect. In addi- 
tion, the calculated lift coefficient during the deep-stall 
phase of the motion is generally too low, and becomes even 
more so as the reduced frequency is increased. It appears 
that a decrease in accuracy results when the values of both 
the reduced frequency and the amplitude of oscillation are 
large . 

5. It is believed that the quantitative accuracy of the 
math model can be substantially improved by incorporating a 
time-delay factor to account for the delay in stall during the 
a > 0 portion of the cycle. Also, based on the observed static 
hysteresis, an improvement in the extent of the dynamic 
hysteresis may be realized by allowing for a separate A func- 
tion during the a < 0 portion of the cycle. It is recom- 
mended that consideration be given to both of these poten- 
tial improvements to the model before examining the drag 
and moment responses. 

Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, California 94035, June 18, 1984 
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Figure 3.- Aeromechanics Laboratory 21- by 31-Centimeter Water Tunnel. 
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Figure 4 - Model installation and balance system for force and moment. 
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Figure 5.- Load contributions to balance measurements. 
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Figure 1 1 Load measurements for a = 10° + 10° sin ojt at k = 0.10. 
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Figure 15.— Continued. 













Figure 15.- Continued. 














Figure 15.— Continued. 
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Figure 15 .— Continued. 










Figure 15.— Continued. 



















Figure 15.— Continued. 







Figure 15.— Continued. 









Figure 15.— Continued. 
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Figure 16.- Concluded. 
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Figure 27.- Concluded. 
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